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A vector lattice model for stresses in granular materials is proposed. A two dimensional pile built 
by pouring from a point is constructed numerically according to this model. Remarkably, the pile 
violates the Mohr Coulomb stability criterion for granular matter, probably because of the inherent 
anisotropy of such poured piles. The numerical results are also compared to the earlier continuum 
FPA model and the (scalar) lattice g-model. 
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Understanding how granular materials sustain exter- 
nally applied loads is an interesting and difficult scien- 
tific problem, with obvious engineering implications. Un- 
like for elastic materials, the slippages that occur when a 
granular heap is built make it impossible to define a dis- 
placement (and thence a strain) field. One thus has to 
solve for a tensor stress field (which cannot be expressed 
in terms of a vector displacement field) using the vector 
force balance equation. In two dimensions, where the 
stress tensor has three independent components, this re- 
sults in one 'missing' equation. In three dimensions, the 
situation is even worse. 

One of the main tools used to deal with this difficulty is 
the Mohr Coulomb (MC) stability criterion. This states 
that the ratio of the shear to the normal stress at any 
point inside a cohesionless granular heap, with any ori- 
entation of axes, cannot exceed the coefficient of friction 
fi. This is natural, since otherwise one would expect an 
appropriately located (and oriented) slip plane to spon- 
taneously destabilize the heap. Using the MC inequality 
and symmetry (or further assumptions), it is sometimes 
possible to obtain useful results on the stability of gran- 
ular heaps Q. For instance, the elastoplastic model ||] 
of granular stresses relies crucially on this criterion. The 
fixed principal axes (FPA) model does not invoke 
the MC criterion, but obtains it as a byproduct. 

Both the elastoplastic and the FPA theory yield the 
surprising result that the MC criterion is saturated in 
a large outer region of a two-dimensional pile built by 
pouring from a point. This implies B that an infinitesi- 
mal perturbation inside the outer region should destabi- 
lize the pile, while one would expect that a poured pile 
should only be unstable at the surface. 

A byproduct of the MC criterion is that the angle of 
repose <f> should satisfy the equation 
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In two dimensions, the MC criterion becomes 
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isotropic. This is a questionable assumption for a granu- 
lar heap built by pouring from a point. As the grains roll 
down the slopes of the heap, anisotropy can be 'frozen 
in', with different behavior parallel and perpendicular to 
the slope. In fact, for a two dimensional heap of identical 
disks in a perfect triangular lattice || , the left hand side 
of Eq. (||) goes to unity at the corner of an infinitely high 
heap, while cf> = 7r/3. The motivation for this paper is to 
investigate whether this anisotropy persists away from 
the special case of the perfect lattice. 




Although the MC criterion seems reasonable, there 
is the implicit assumption that the granular medium is 



(a) 



FIG. 1. (a) Triangular lattice. Each site rests on its two 
neighbors below, and receives forces from its neighbors above, 
(b) The incoming force is too horizontal, and the site is unsta- 
ble. The downwards bond going left is broken, and a horizon- 
tal bond is established. The angles 0i and 62 are independent 
random variables at each site. 

In this paper, we consider a two dimensional lattice 
model, with the grains placed on a triangular lattice (see 
Figure Each grain rests on its two neighbors in the 
row below, and supports its neighbors in the row above. 
Adjoining sites in the same horizontal row are not con- 
nected. Forces propagate downwards from each grain to 
its descendants. In addition to the force from its prede- 
cessors, a grain also transmits its own weight downwards. 
This weight is taken to be the same for every grain. 

In order to include disorder in the model, the outgoing 
bond angles 9\ and #2 from a site are chosen randomly. 
For N lattice sites, there are 2N random angles, each of 
which is chosen independently from a uniform distribu- 
tion over the interval [# m i n , max ] ■ The width of this in- 
terval is a measure of the degree of disorder. With equal 
sized grains, even if the lattice is viewed as a topologi- 
cal rather than a geometrical structure, the bond angles 
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would not be all independent, but such correlations are 
ignored for simplicity. 

Vector force balance is enforced on every lattice site. 
Solving for the inter-grain force network, proceeding 
downwards from row to row, inevitably leads to occa- 
sional negative (i.e. attractive) inter-grain forces. This 
is unphysical for a cohesionless granular medium. Ac- 
cordingly, if the force between a grain and one of its de- 
scendants is negative, we interpret this to mean that the 
grain is unstable, and rolls in the direction of the other 
descendant. Within the model, this is accommodated by 
breaking the bond with negative force, and establishing a 
new bond with the adjacent site (in the same row as the 
unstable grain) in the opposite direction (see Figure [l]) . 
The new bond established necessitates recomputing force 
balance for the sites in the row; several iterations can be 
required before a stable configuration is achieved. 

In reality, a moved grain would also change the an- 
gles for its unbroken bonds. Moreover, since the pile is 
not really a regular lattice, the link established to the 
adjacent site would not be horizontal. Both of these is- 
sues are ignored. The idea of unstable grains moving and 
readjusting their contacts has been suggested earlier |Q, 
although the implementation of this process and the spe- 
cific system studied are different here. 

The model described above is similar to the earlier 
g-model with one key difference: the q-model has 
scalar forces, only keeping track of the vertical force on 
each grain. Even with this limitation, the g-model suc- 
cessfully accounts for the experimental result P, p!(i| , pT| 
that the vertical forces on individual grains at the bot- 
tom of a granular heap have a distribution P(f) that 
decays exponentially for large /. It also agrees with the 
experimental observation that the vertical forces fa 
and fj on two grains at the bottom of a heap satisfy 
(fifj) — ifi)(fj)- However, the absence of vector forces 
in the g-model prevents it from reproducing the visually 
most striking feature of granular forces: the existence of 
"force chains" Jl2|,|| , a sparse network of grains that ex- 
perience large stresses. Further, while forces on a single 
grain level are important for the failure of granular ma- 
terials, one is often interested in a more coarse grained 
continuum description. There have been earlier work 
on vectorizing the q- model indirectly [ [""""jf , and directly 
with a method different from the one here H] . We shall 
see that the vector lattice model proposed in this paper 
agrees well with the same experimental results as the q- 
model, and in addition yields force chains, corroborates 
some aspects of the FPA model while disagreeing with 
others, and permits a test of the MC stability criterion. 

The specific system we consider using the vector lat- 
tice model is that of a granular pile built by pouring 
from a point. This is because this is the system on which 
most of the FPA model analysis has concentrated, and 
also because we have reasons to question the validity of 
the MC criterion in such a system. For convenience, in 



the numerical simulations, the pile is built from the top 
down instead of from the bottom up. Thus one starts 
with one lattice site in the top row and computes the 
outgoing forces, proceeding to two sites in the next row 
and so on. If a row has m sites, the next row has m + 1 
sites unless the sites at the edge require adjacent sites 
to establish horizontal bonds with. All the numerics in 
this paper are with bond angles ranging over [71-/6, 7r/3], 
unless specified otherwise. 




FIG. 2. A pile 400 rows deep. The outer lines are the 
boundaries of the pile, and show that it has a well defined 
slope. In any row, bonds with a force greater than 4 times 
the average bond force in that row are marked, and show di- 
rected force chains. The thin straight lines are the boundaries 
of the inner region of the pile. The force chains are concen- 
trated in the inner region. (If the threshold were less than 4, 
there would be more chains in the outer region.) 

Figure [""] shows the result of such a simulation, with 
400 rows. One first notices that there is a well defined 
slope to the pile. From the manner in which the pile is 
built up, we view the central m sites in the m'th row as 
the interior of the pile, and the sites that flank this inner 
region as buttresses, necessary to stabilize the pile. This 
division into an inner and outer region is in accordance 
with the FPA model. With the bonds that support un- 
usually high stresses marked, one also sees force chains 
concentrated in the inner region of the pile, in accor- 
dance with experiment ]f5| . Although the force chains 
seem perfectly straight, this is because the sites are on 
a regular lattice. With a more realistic model, where 
the random bond angles would be accompanied by an ir- 
regular lattice, the chains should meander, with a mean 
orientation of (6). 

If the horizontal and vertical forces across the bottom 
of the pile (averaged over many runs) are plotted, the 
vertical force f y has a flat peak in the inner region and 
falls off steeply in the outer region. The variation in f y 
across the inner region is less than ±5%. The horizon- 
tal force f x varies linearly across the inner region, peaks 



2 



roughly at the boundary between the inner and outer re- 
gions, and falls off in the outer region. The linear profile 
for f x shows that the pile has a strong tendency to splay. 
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FIG. 3. (a) Components of the stress tensor across the 
bottom of the pile. The boundary of the inner region is 
at x = ±4000. (b) The ratio a xx /a yy from the same data, 
showing that a yy is not proportional to a xx . (There is a slight 
sharpening of the peaks with increasing pile size.) a xx cannot 
be expressed even as a linear combination of o yy and |(Xxy|. 



One can coarse grain and compute 'continuum' 
stresses: Figure || shows the components of the stress 
tensor, calculated by coarse graining over twenty lattice 
sites at the bottom of a pile of 8000 rows, and then av- 
eraging over 2000 runs. Not surprisingly, the division 
between the inner and outer regions seen in the force 
components f x and f y is also seen here. 

The slight curvature seen in a yy in the central region 
does not decrease as a function of the size of the pile. This 
is accentuated in the plot of a xx /a yy . The steady increase 
in this stress ratio as one moves away from the symmetry 
axis violates the BCC hypothesis || of a yy oc a xx ; this 
cannot be cured by using the FPA or OSL hypotheses Q 
of <J XX = k\(jyy + k2\<J xy \, since a xx ja yy is quadratic at 
the center of the pile, whereas \o xy \ has a linear kink. 
More importantly, the violation of the FPA hypothesis 
can be understood physically. At least in the inner re- 
gion, the splay in the pile makes it increasingly likely 



for sites to be destabilized outwards as one moves away 
from the axis of symmetry. Such destabilized sites es- 
tablish bonds with their adjacent sites. Thus while the 
forces emerging from a site in the central region prop- 
agate downwards to its descendants, there is a steady 
drift towards more horizontal force propagation as one 
moves outwards. This reorients the principal axes of the 
stress tensor and increases a xx /a yy . Changing the range 
of bond angles, [# m in, #max] has no qualitative effect on 
the results discussed so far: the size of the buttresses 
flanking the inner region naturally increases when the 
bond angle range is increased. 
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FIG. 4. Semi-log plot of P(f) vs / from the inner region 
of piles of varying depths; / is in arbitrary units and scaled 
with system size. The asymptote is roughly straight, and 
grows straighter with increasing system size. The inset (axis 
label on right) shows a similar plot with bond angles almost 
spanning [0, n/2]. 

The approximately flat profile of the vertical force in 
the inner region allows us to compute the distribution 
of vertical forces in this region. Figure |^ was obtained 
by averaging over many runs and lattice sites, but the 
forces on a single site in a single run have large fluctua- 
tions. Figure |5] shows a semi-log plot of the probability 
P(f) for a grain in (the inner region of) the bottom row 
to experience a force /. Although there is a slight devi- 
ation from an exponential for large /, this steadily de- 
creases as the size of the pile is increased, suggesting an 
exponential tail to P(f) for infinite system size. If the 
range of bond angles is reduced, the evolution towards 
an exponential tail for P(f) with increasing system size 
is slower, but is nevertheless still manifest. Conversely, 
the inset to Figure [| shows lnP(/) with a bond angle 
range of [0.0l7r/3, 1.497r/3]. Even though the pile is only 
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400 rows deep, P(f) clearly has an exponential tail. 
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FIG. 5. Mohr Coulomb ratio across the bottom of the pile. 
The components of the stress tensor are calculated as in Fig. ^| 
The dashed line is sin 2 a, where a is the angle of repose. The 
curve clearly goes above the line, showing that Eq. (Q) is vio- 
lated. There is a slight rounding of the peaks with increasing 
pile depth. 



sin 2 a, where a is the angle of repose. It is slightly tricky 
to define the angle of repose when a regular lattice is 
combined with random bond angles. If the width of the 
pile at a depth of m rows is w(m), and the large m limit 
of m/w(m) is A, since the average bond angle is 7r/4, a 
is taken to be tan -1 A. The numerical curve is seen to 
clearly go above the straight line, denoting a violation of 
the MC criterion. 

In this paper, we have analyzed a new two dimen- 
sional lattice model with vector forces for granular ma- 
terials. On an individual grain level, the model yields 
force chains, arching, and is consistent with an exponen- 
tial tail to the distribution of vertical forces. Continuum 
stresses show an inner and outer region for a pile poured 
from a point, but the principal axes of the stress tensor 
are not fixed. Most importantly, the model violates the 
Mohr Coulomb criterion, calling in question the validity 
of this criterion for poured granular heaps. 

I thank Sid Nagel, Deepak Dhar and Narayanan Menon 
for useful discussions. 



The small / behavior of P{f) depends on the range 
of bond angles. However, in all cases one sees that 
P(f — > 0) 7^ 0. This is in accordance with experi- 
ments [pd| , but not with the (/-model. In the model in 
this paper, if a lattice site receives a large force from 
above, it transmits large forces to both its descendants. 
Since it is unlikely that either of its descendants gets a 
large force from its other ancestor, both descendants are 
usually destabilized outwards. The common descendant 
of the two destabilized grains thus receives no force from 
either of its ancestors. This yields a <5-function in P(f) 
at / = (more precisely, at / equal to the weight of a 
single grain), which leads to P(f —> 0) ^ 0. Thus the vec- 
tor lattice model agrees with experiments in this respect 
because it has arching, which is important in granular 
materials but is absent in the g-modcl. 

One can also compute the correlation between vertical 
forces on adjacent sites (in the inner region of the pile) 
in the same horizontal row Unlike the g-model this is not 
exactly zero, and depends on the bond angle distribu- 
tion, but is typically about a few percent of the variance 
in the vertical forces. This is reasonably consistent with 
the experimental observation jll]] that the vertical forces 
on different grains are uncorrelated. 

Finally, Figure || shows the MC ratio, the left hand 
side of Eq.(0), obtained from the coarse-grained stresses 
of Figure ||T The dashed line shown is at a height of 



[1] R.M. Nedderman, Statics and Kinematics of Granular 

Materials (Cambridge Univ. Press, Cambridge, 1992). 
[2] F. Cantelaube and J.D. Goddard, in Powders and Grains 

97, R.P. Behringer and J.T. Jenkins eds. (Balkema 1997). 
[3] J.-P. Bouchaud, M.E. Cates and P.J. Claudin, J. Phys. I 

France 5, 639 (1995). 
[4] J. P. Wittmer, M.E. Cates and P. Claudin, J. Phys. I 

France 7, 39 (1997). 
[5] O. Narayan and S.R. Nagel, Physica A 264, 75 (1999). 
[6] D.C. Hong, Phys. Rev. E 47, 760 (1993). 
[7] A.V. Tkachenko and T.A. Witten, unpublished. 
[8] C.-H. Liu, S.R. Nagel, DA. Schecter, S.N. Coppersmith, 

S.N. Majumdar, O. Narayan and T.A. Witten, Science 

269, 513 (1995). 
[9] S.N. Coppersmith, C.-H. Liu, S. Majumdar, O. Narayan 

and T.A. Witten, Phys. Rev. E53, 4673 (1996). 
[10] B. Miller, C. O'Hern and R.P. Behringer, Phys. Rev. 

Lett. 77, 3110 (1996). 
[11] D.M. Mueth, H.M. Jaeger and S.R. Nagel, Phys. Rev. 

E57, 3164 (1998). 
[12] P. Dantu, Ann. Ponts Chaussees IV, 144 (1967); A. 

Drescher and G. de Josselin de Jong, J. Mech. Phys. 

Solids 20, 337 (1972); T. Travers, M. Ammi, D. Bideau, 

A. Gervois, J.C. Messager and J. P. Troadec, Europhys. 

Lett. 4, 329 (1987). 
[13] P. Claudin and J-P. Bouchaud, Phys. Rev. Lett. 78, 231 

(1997). 

[14] J.E.S. Socolar, Phys. Rev. E57, 3204 (1998). 

[15] D.W. Howell, J. Ceng and R.P. Behringer, unpublished. 



4 



